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Abstract 

Parallel MRI is a fast imaging technique that enables the acquisition of highly resolved images in space or/and 
in time. The performance of parallel imaging strongly depends on the reconstruction algorithm, which can proceed 
either in the original /c-space (GRAPPA, SMASH) or in the image domain (SENSE-like methods). To improve 
the performance of the widely used SENSE algorithm, 2D- or slice-specific regularization in the wavelet domain 
has been deeply investigated. In this paper, we extend this approach using 3D-wavelet representations in order 
to handle all slices together and address reconstruction artifacts which propagate across adjacent slices. The 
gain induced by such extension (3D-Unconstrained Wavelet Regularized -SENSE: 3D-UWR-SENSE) is validated 
on anatomical image reconstruction where no temporal acquisition is considered. Another important extension 
accounts for temporal correlations that exist between successive scans in functional MRI (fMRI). In addition to the 
case of 2D+t acquisition schemes addressed by some other methods like At-FOCUSS, our approach allows us to 
deal with 3D+t acquisition schemes which are widely used in neuroimaging. The resulting 3D-UWR-SENSE and 
4D-UWR-SENSE reconstruction schemes are fully unsupervised in the sense that all regularization parameters are 
estimated in the maximum likelihood sense on a reference scan. The gain induced by such extensions is illustrated 
on both anatomical and functional image reconstruction, and also measured in terms of statistical sensitivity 
for the 4D-UWR-SENSE approach during a fast event-related fMRI protocol. Our 4D-UWR-SENSE algorithm 
outperforms the SENSE reconstruction at the subject and group levels (15 subjects) for different contrasts of 
interest (eg, motor or computation tasks) and using different parallel acceleration factors (R = 2 and R = 4) on 
2 x 2 x 3mm 3 EPI images. 



1 Introduction 

Reducing scanning time in Magnetic Resonance Imaging (MRI) exams remains a worldwide chal- 
lenging issue. The expected benefits of a faster acquisition can be summarized as follows: i.) limit 
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patient's exposure to the MRI environment either for safety or discomfort reasons, ii.) maintain a 
strong robustness in the acquisition with respect to subject's motion artifacts, Hi.) limit geometric 
distortions or maintain high image quality and iv) acquire more spatially or temporally resolved 
images in the same or even reduced amount of time [1,2]. The basic idea to make MRI acquisitions 
faster or to improve spatial resolution at fixed scanning time consists of reducing the amount of 
acquired samples in the fc-space and developing dedicated reconstruction pipelines. To achieve this 
goal, three main research avenues have been developed so far: 

• parallel imaging that relies on a geometrical principle involving multiple receiver coils with 
complementary sensitivity profiles. This enables the reduction of the number of fc-space lines 
to be acquired without degrading spatial resolution or truncating the Field-Of-View (FOV), 
and thus requires the unfolding of reduced FOV coil-specific images to reconstruct the full 
FOV image [3-5]. 

• Compressed Sensing (CS) MRI that exploits the implicit sparsity in MR images to signifi- 
cantly undersample the fc-space and randomly select incoherent samples regarding their spec- 
tral contribution to the MR image [6]. This approach remains usable with birdcage coil and 
does not require parallel imaging while it can be combined with this acquisition strategy [7]. 

• In the context of dynamic MRI, fast parallel acquisition schemes have been proposed to 
increase the acquisition rate by reducing the amount of acquired fc-space samples in each frame 
using interleaved partial fc-space sampling between successive frames (UNFOLD approach [8]). 
To further reduce the scanning time, a strategy named fct-BLAST taking advantage of both 
the spatial (actually in the fc-space) and temporal correlations between successive scans in 
the dataset has been pushed forward [9]. 

In parallel Magnetic Resonance Imaging (pMRI), many reconstruction methods like SMASH (Si- 
multaneous Acquisition of Spatial Harmonics) [3], GRAPPA (Generalized Autocalibrating Partially 
Parallel Acquisitions) [5] and SENSE (Sensitivity Encoding) [4] have been proposed in the literature 
to reconstruct a full FOV image from multiple fc-space undersampled images acquired on separate 
channels. The main difference between these two classes of methods lies in the space on which they 
operate. GRAPPA performs multichannels full FOV reconstruction in the fc-space domain whereas 
SENSE carries out the unfolding process in the image domain: all the undersampled images are 
first reconstructed by inverse Fourier transform before combining them to unwrap the full FOV 
image. Another difference is that GRAPPA is autocalibrated, while SENSE needs a separate coil 
sensitivity estimation step based on a reference scan. Note however that an autocalibrated version 
of SENSE is also available for instance in Siemens scanners and called mSENSE hereafter. SENSE as 
well as GRAPPA methods may suffer from strong artifacts when high values of acceleration factor 
R are considered in the imaging setup or when they are applied to Echo Planar Imaging (EPI) , the 
sequence involved in fMRI experiments. 

In the dynamic MRI context, combined strategies mixing parallel imaging or pMRI and accel- 
erated sampling schemes along the temporal axis have also been investigated. The corresponding 
reconstruction algorithms have been referenced to as H-SENSE [9,10], H-GRAPPA [11]. More 
recently, optimized versions of H-BLAST and H-SENSE reconstruction algorithms referenced to 
as H-FOCUSS [12, 13] have been designed to combine the CS theory in space with Fourier or alter- 
native transforms along the time axis. They enable to further reduce data acquisition time without 
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significantly compromising image quality, provided that the image sequence exhibits a high degree 
of spatio-temporal correlation, either by nature or by design. Typical examples that enter in this 
context are i.) dynamic MRI capturing an organ (liver, kidney, heart) during a quasi-periodic 
motion due to the respiratory cycle and cardiac beat and ii.) functional MRI based on periodic 
blocked design. 

Compared to their mSENSE-like alternative where the centre of the fc-space is acquired only for 
the first repetition, the limitation of these methods remains their necessity to acquire the central 
fc-space area at each repetition along the time axis, which decreases the acceleration factor. Also, 
these methods remain inapplicable for 3D acquisition sequences where a 3D encoding scheme is 
used to acquire the whole brain volume. Such methods cannot also be used when only static MRI 
acquisitions (for example anatomical image) are considered by turning off the temporal repetition. 

Moreover, the interleaved partial fc-space sampling cannot be exploited in aperiodic dynamic 
acquisition schemes, even if a certain amount of serial correlation may be exhibited between suc- 
cessive scans. Typical aperiodic situations arise in resting state or paradigm-free fMRI as well as 
during fast-event related fMRI paradigms [14,15]. In the first case, ongoing spontaneous brain 
activity is recorded without any experimental design in order to probe functional connectivity or 
resting-state networks [14, 16, 17]. In the second case, the presence of jittering in the experimen- 
tal design introduces a trial- varying delay between the stimulus and acquisition time points [18]. 
Moreover, the use of several stimulus types presented in a random order prevents the use of an 
interleaved fc-space sampling strategy along the temporal axis. Indeed, in such circumstances there 
is no guarantee that the BOLD response is quasi-periodic. In addition, the vast majority of fMRI 
studies in neurosciences make use either of resting state fMRI or fast event-related designs [18, 19]. 
In these two contexts, the most reliable acquisition strategy is the "scan and repeat" approach, 
while it is suboptimal. To the best of our knowledge only one kt- contribution (fct-GRAPPA [11]) 
has claimed its ability to accurately reconstruct fMRI images in aperiodic paradigms. 

Overview of our contribution 

The present paper therefore aims at proposing a new 4-dimensional parallel imaging reconstruc- 
tion algorithm that can be adopted irrespective of the nature of the encoding scheme or the fMRI 
paradigm. In particular, we show here that our approach outperforms classical algorithms like 
SENSE for anatomical images reconstruction and also when using a fast event-related design in- 
volving multiple stimulus types in a random and thus aperiodic order. 

In the fMRI literature, a few studies have been conducted to measure the impact of the par- 
allel imaging reconstruction algorithm on EPI volumes and subsequent statistical sensitivity for 
detecting evoked brain activity [2,20-22]. In these works, reliable activations were detected for 
an acceleration factor up to 3. More recently, a special attention has been paid in [23] to assess 
the performance of dynamic MRI reconstruction algorithms on BOLD fMRI sensitivity. These au- 
thors have reported that fct-based approaches perform better than conventional SENSE for BOLD 
fMRI in the sense that reliable sensitivity may be achieved at higher undersampling factors (up to 
5). However, most of the time, these comparisons are made on a small group of individuals and 
statistical analysis is only performed at the subject level. Here, we perform the comparison of sev- 
eral pMRI reconstruction algorithms both at the subject and group levels for different acceleration 
factors. 

In neuroimaging, reconstruction artifacts can drastically disturb subsequent statistical analysis 
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such as brain activation detection. Regularized SENSE methods have been proposed in the liter- 
ature to improve the robustness of the solution [24-28]. Some of them apply quadratic or Total 
Variation (TV) regularizations while others resort to 2D regularization in the wavelet transform 
domain (e.g. UWR-SENSE [29]). The latter strategy has proved its efficiency on the reconstruction 
of anatomical or functional (resting-state only) data, compared to standard SENSE and TV-based 
regularization [27,29]. More recently, UWR-SENSE has been assessed on EPI images and com- 
pared with mSENSE on the same data acquired during a brain activation fMRI experiment [30]. This 
comparison was performed at the subject level on a few subjects only. 

Besides, except some non-regularized contributions like 3D GRAPPA [31], most of the avail- 
able reconstruction methods in the literature operate slice by slice and thus reconstruct each slice 
irrespective of its neighbours. Iterating over slices is thus necessary to get the whole 3D volume. 
This observation led us to consider 3D or whole brain reconstruction as a single step in which all 
slices are treated together by making use of 3D wavelet transforms and 3D sparsity promoting 
regularization terms in the wavelet domain. Following the same principle, an fMRI run usually 
consists of several tens of successive scans that are reconstructed independently one to another and 
thus implicitly assumed to be independent of each other. Iterating over all acquired 3D volumes 
remains the classical approach to reconstruct the 4D or 3D + t dataset associated with a fMRI run. 
However, it has been shown for a long while that fMRI data are serially correlated in time even 
under the null hypothesis (i.e., ongoing activity only) [32-34]. To capture this dependence between 
successive time points, an autoregressive model has demonstrated its relevance [35-38], especially 
when its parameters and optimal order can vary in space for instance with tissue type [39,40]. 
Hence, it makes sense to account for this temporal structure at the reconstruction step. 

These two key ideas play a central role in the present paper that extends the UWR-SENSE 
approach [29] through a more general convex objective function to be minimized that accounts for 
3D spatial and temporal dependencies between successive slices and acquisition times, respectively. 
Our novel criterion relies on a fidelity data term that combines all time points or repetitions and 
involves a 3D (and not only 2D as in [29]) wavelet transform. Associated with regularization in the 
wavelet domain, an additional regularization term along the temporal dimension of the 4D dataset 
enables to account for the temporal correlation in the time series. Compared to other methods like 
H-FOCUSS which can only deal with 2D+t acquisition (slice by slice), our SENSE-based extension 
allows us to deal with the whole brain volume at once while keeping higher acceleration factors 
since the fc-space centre is not acquired at every repetition. 

The development of the proposed method (named 4D-UWR-SENSE) was made possible due to 
recent advances in nonsmooth convex optimization. Indeed, it is based on the Parallel ProXimal 
Algorithm (PPXA) which can address a broader scope of optimization problems than the forward- 
backward and Douglas-Rachford methods employed in [29] , or even Fast Iterative Soft Thresholding 
Algorithm (FISTA)-like as used in [41]. All these algorithms are only able to optimize the sum of 
two convex functions, while the PPXA allows the optimization of any sum (even more than two) 
of convex functions. Indeed, using such an algorithm is necessary in our case since as detailed 
in Section 4, our function to be optimized is made up of more than two convex functions due 
to the additional temporal regularization. For this reason, the proposed reconstruction method 
cannot rely on more standard optimization algorithms like in [29, 41]. Moreover, our work can also 
be viewed as an extension of the approach in [41] to a fully 4D framework. In the latter paper, a 
preconditioned version of FISTA was proposed in order to minimize a spatially regularized criterion 
involving the t\ norm of the wavelet coefficients of the target image. 
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4D-UWR-SENSE leads to improvements in the retrieval of a reliable Signal-to-Noise Ratio (SNR) 
between the acquired volumes and also in the enhancement of the detection of BOLD effects or 
evoked activations that occur in response to the delivered stimuli during the fMRI experiment. The 
present paper therefore aims at demonstrating that the 4D-UWR-SENSE approach outperfoms its 
SENSE-like alternatives not only in terms of reduced artifacts, but also in terms of statistical 
sensitivity at the voxel and cluster levels, in intra-subject and group studies. 

The rest of this paper is organized as follows. Section 2 recalls the general parallel MRI framework. 
We then describe the proposed reconstruction algorithm in Section 3 before illustrating related 
experimental results in Section 5. Finally, some discussions and conclusions are drawn in Section 6. 

2 Parallel imaging in MRI 

In parallel MRI, an array of L coils is employed to measure the spin density ~p into the object 
under investigation 1 . The signal di received by each coil £ (1 < £ < L) is the Fourier transform 
of the desired 2D field ~p E M Xxy on the specified FOV weighted by the coil sensitivity profile 
evaluated at some location k r = (k x , k y ) T in the fc-space: 



where ni(k r ) is a coil-dependent additive zero-mean Gaussian noise, which is independent and 
identically distributed (iid) in the fc-space, and r = (x, y) T E X x Y is the spatial position in the 
image domain (- T being the transpose operator). The size of the reduced FOV acquired data dn 
in the fc-space clearly depends on the sampling scheme. In contrast with cardiac imaging where 
significant heart motion makes suitable the use of spiral acquisition schemes, a Cartesian coordinate 
system is generally adopted in the neuroimaging context. In parallel MRI, the sampling period 
along the phase encoding direction is R times larger than the one used for conventional acquisition, 
R < L being the reduction factor. To recover full FOV images, many algorithms have been proposed 
but only SENSE-like [42] and GRAPPA-like [5] methods are provided by scanner manufacturers. 
In what follows, we focus on SENSE-like methods operating in the image domain. 

Let Ay = be the aliasing period and y the position in the image domain along the phase 
encoding direction. Let x be the position in the image domain along the frequency encoding 
direction. A 2D inverse Fourier transform allows us to recover the measured signal in the image 
domain. By accounting for the fc-space undersampling at iterate, the inverse Fourier transform 
gives us the spatial counterpart of Eq. (1) in matrix form: 




(i) 



d(r) = S(r)p(r) + n(r), 



(2) 



where 




si(x,y) 



Sl (x,y + (R-l)Ay) 



n(r) — 



ni{x,y) 
n 2 {x,y) 



SL{x,y) 



s L (x,y + (R-l)Ay) 



nL{x,y) 



lr The overbar is used to distinguish the "true" data from a generic variable. 
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p(r) = 



p(x, y + Ay) 
p(x,y + (R-l)Ay) 



A 



and d(r) — 



di(x,y) 
d 2 (x,y) 

d L (x,y) 



(3) 



Based upon this model, the reconstruction step consists of solving Eq. (2) so as to recover p(r) 
from d(r) and an estimate of S(r) at each spatial position r = (x,y) J . The spatial mixture or 
sensitivity matrix S(r) is estimated using a reference scan and varies according to the coil geometry. 
Note that the coil images (^)kkl as well as the sought image ~p are complex- valued, although |p| 
is only considered for visualization purposes. The next section describes the widely used SENSE 
algorithm as well as its regularized extensions. 

3 Reconstruction algorithms 
3.1 ID-SENSE 

In its simplest form, SENSE imaging amounts to solving a one-dimensional inversion problem due 
to the separability of the Fourier transform. Note however that this inverse problem admits a 
two-dimensional extension in 3D imaging sequences like Echo Volume Imaging (EVI) [2] where un- 
dersampling occurs in two fc-space directions. The ID-SENSE reconstruction method [42] actually 
minimizes a Weighted Least Squares (WLS) criterion Jwls given by: 



Jwls(p) 



E 

re{l,...,X}x{l, 



d(r) - S(r)p(r) 



|2 



(4) 



,Y/R} 



where 



I*- 



1 = J(-) H * -1 (-)> and the noise covariance matrix \l/ is usually estimated based on 
L acquired images (d^)i<£<L from all coils without radio frequency pulse. Hence, the SENSE full 
FOV image is nothing but the maximum likelihood estimate under Gaussian noise assumption, 
which admits the following closed- form expression at each spatial position r: 



PwlsW = (S^r^Sir^S^r^dir), 



(5) 



where (-) H (respectively (•)") stands for the transposed complex conjugate (respectively pseudo- 
inverse). It should be noticed here that the described ID-SENSE reconstruction method has been 
designed to reconstruct one slice (2D image). To reconstruct a full volume, the ID-SENSE recon- 
struction algorithm has to be iterated over all slices. 

In practice, the performance of the SENSE method is limited because of i) different sources of 
noise such as distortions in the measurements d(r) and the estimation of S(r) mainly at brain/air 
interfaces, and ii) the ill-conditioning of S(r). To enhance the robustness of the solution to this 
ill-posed problem, a regularization is usually introduced in the reconstruction process. To improve 
results obtained with quadratic regularization techniques [24,25], edge-preserving regularization 
has been widely investigated in the pMRI reconstruction literature. For instance, reconstruction 
methods based on Total Variation (TV) regularization have been proposed in a number of recent 
works like [43,44]. However, TV is mostly adapted to piecewise constant images, which are not 
always accurate models in MRI, especially in fMRI. As investigated by Chaari et al. [29], Liu et 
al. [28] and Guerquin-Kern et al. [41], regularization in the Wavelet Transform (WT) domain is a 
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powerful tool to improve SENSE reconstruction. In what follows, we summarize the principles of 
the wavelet-based regularization approach. 

3.2 Proposed wavelet regularized SENSE 

Akin to [29] where a regularized reconstruction algorithm relying on 2D separable WTs was in- 
vestigated, to the best of our knowledge, all the existing approaches in the pMRI regularization 
literature proceed slice by slice. The drawback of this strategy is that no spatial continuity between 
adjacent slices is taken into account since the slices are processed independently. Moreover, since 
the whole brain volume has to be acquired several times in an fMRI study, iterating over all the 
acquired 3D volumes is then necessary in order to reconstruct a 4D data volume corresponding to 
a fMRI session. 

Consequently, the 3D volumes are supposed independent whereas fMRI time-series are serially 
correlated in time because of two distinct effects: the BOLD signal itself is a low-pass filtered version 
of the neural activity, and physiological artifacts make the fMRI time series strongly dependent. 
For these reasons, modeling temporal dependence across scans at the reconstruction step may 
impact subsequent statistical analysis. This has motivated the extension of the wavelet regularized 
reconstruction approach in [29] in order to: 

• account for 3D spatial dependencies between adjacent slices by using 3D WTs, 

• exploit the temporal dependency between acquired 3D volumes by applying an additional 
regularization term along the temporal dimension of the 4D dataset. 

This additional regularization will help in increasing the Signal to Noise Ratio (SNR) through 
the acquired volumes, and therefore enhance the reliability of the statistical analysis in fMRI. These 
temporal dependencies have also been used in the dynamic MRI literature in order to improve 
the reconstruction quality in conventional MRI [45]. However, since the imaged object geometry 
in the latter context generally changes during the acquisition, taking into account the temporal 
regularization in the reconstruction process is very difficult. 

To deal with a 4D reconstruction of the N r acquired volumes, we will first rewrite the observation 
model in Eq. (2) as follows: 

d t {r) = S{r)p\r)+n t {r), (6) 

where t E {l,...,7V r } is the acquisition time and r — (x,y,z) is the 3D spatial position, z E 
{1, . . . , Z} being the position along the third direction (slice selection one). 

At a given time t, the full FOV 3D complex- valued image ~p l of size X xY x Z can be seen as an 
element of the Euclidean space C K with K = X xY x Z endowed with the standard inner product 
(• | •) and norm || • || . We employ a dyadic 3D orthonormal wavelet decomposition operator T over j max 
resolution levels. The coefficient field resulting from the wavelet decomposition of a target image p l 
is defined as <? = (C*„, (<o,i)oeO,i<i< Jmax ) with o £ O = {0, 1} 3 \{(0, 0, 0)}, £ = (£ >k )i<k<K jinax and 
Co j — (Co j k)t<k<Kj where Kj — K2~ 3 i is the number of wavelet coefficients in a given subband at 
resolution j (by assuming that X, Y and Z are multiple of 2 Jmax ). Adopting such a notation, the 
wavelet coefficients have been reindexed so that Ca denotes the approximation coefficient vector at 
the resolution level j m ax 5 while Coj denotes the detail coefficient vector at the orientation o and 
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resolution level j. Using 3D dyadic WTs allows us to smooth reconstruction artifacts along the 
slice selection direction, which is not possible using a slice by slice operating approach. 

The proposed regularization procedure relies on the introduction of two penalty terms. The 
first penalty term describes the prior 3D spatial knowledge about the wavelet coefficients of the 
target solution and it is expressed as: 

5(0 = E [ E + E E E *oA£j,k)\ , (7) 

t=l k=l oeO j=l k=l 

where ( = (C 1 , £ 2 , . . . , ( Nr ) and we have, for every o E O and j E {1, . . . , jmax} (and similarly for 
<3> a relative for the approximation coefficients), 

V£eC, * o>J -(0 = *S(0 + *S(0 (8) 

where <&^(0 = < e |Re(£ - Mo , ;)| + %*|Ite(£-/x j)| 2 and ^(C) = <|Im(£ - Mo j)| + % i |Im(^- 
Mo,j)| 2 with = /i^ + z/i*™ E C, and a^, /3^|, a£™, /J*™ are some positive real constants. 
Hereabove, Re(-) and Im(-) (or - Re and - Im ) stand for the real and imaginary parts, respectively. 
For both real and imaginary parts, this regularization term allows keeping a compromise between 
sparsity and smoothness of the wavelet coefficients due to the t\ and £2 terms, respectively. 
The second regularization term penalizes the temporal variation between successive 3D volumes: 

N r 

Ho = Kj2\\ T *?- T ^ t ~ 1 \\ p p ( 9 ) 

t=2 

where T* is the 3D wavelet reconstruction operator. The prior param- 

eters a oJ = (ag,o£}), oJ = H oJ = (f%,V%), « G 

[0, +oo[ and p E [l,+oo[ are unknown and they need to be estimated 
(see Appendix). The used £ p norm gives more flexibility to the temporal penalization term 
by allowing it to promote different levels of sparsity depending on the value of p. Such a penaliza- 
tion has been chosen based on empirical studies that have been conducted on the time-course of 
the BOLD signal at the voxel level. 

The operator T* is then applied to each component ( f of ( to obtain the reconstructed 3D 
volume p t related to the acquisition time t. It should be noticed here that other choices for the 
penalty functions are also possible provided that the convexity of the resulting optimality criterion 
is ensured. This condition enables the use of fast and efficient convex optimization algorithms. 
Adopting this formulation, the minimization process plays a prominent role in the reconstruction 
process, as detailed in Section 4. 



4 Optimization procedure for the 4D reconstruction 

Based on the formulation hereabove, the criterion to be minimized can be written as follows: 

Jst(C) = ^twls(C) + 9(0 + h(0 (10) 
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where Jtwls is defined as 

N r N r 

Jtwls(C) = E^wls(C*) = E E ||d t (r)-5(r)(2*C')(r)||i-i. (11) 

t=l t=l r6{l,...,X}x{l,...,y/fl}x{l,...,Z} 

The minimization of J7st * s performed by resorting to the concept of proximity operators [46], 
which was found to be fruitful in a number of recent works in convex optimization [47-49]. In what 
follows, we recall the definition of a proximity operator: 

Definition 4.1 [46] Let To(x) be the class of lower semicontinuous convex functions from a separa- 
ble real Hilbert space x to ] — oo, +oc] and let tp E Tq(x)- F° r every xGx, the function tp+ || • — x|| 2 /2 
achieves its infimum at a unique point denoted by prox^x. The operator prox^ : x ~~ ^ X ^ s the 
proximity operator of \p. 

In this work, as the observed data are complex- valued, the definition of proximity operators is 
extended to a class of convex functions defined for complex- valued variables. For the function 

3>: C K ^} -oc,+oc] (12) 
x ^ Re (Re(x)) + Im (Im(x)), 

where Re and Im are functions in ro(M x ) and Re(x) (respectively Im(x)) is the vector of the 
real parts (respectively imaginary parts) of the components of x E C x , the proximity operator is 
defined as 

prox^: C K ^C K (13) 
x H> prox^Re(Re(x)) + zprox^i m (Im(x)). 

Let us now provide the expressions of proximity operators involved in our reconstruction prob- 
lem. 



4.1 Proximity operator of the data fidelity term 

According to standard rules on the calculation of proximity operators [49, Table 1.1], the prox- 
imity operator of the data fidelity term Jwls is given for every vector of coefficients ( f (with 
t E {1, ...,JV r }) by proxj = IV, where the image u* is such that Vr E {1,...,X} x 

{i,...,y/i?}x{i,...,z}, 

u\v) = (/ i , + 2S H (r)*- 1 S(r))" 1 (p t (r) + 2S H (r)*- 1 d t (r)), (14) 

where p l = T% f . 
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4.2 Proximity operator of the spatial regularization function 

According to [29], for every resolution level j and orientation o, the proximity operator of the 
spatial regularization function <fr OJ is given by 

V£ € C, prox* o / = si ^^-^)) max{)Re(c _ ^ . } | _ a Re > Q} 

sign(Im(f - /z j)) Im 
+ Z flim _i_ i max{|Im(C - fi ,j)\ ~ a j, 0} + Mo,j (15) 



where the sign function is defined as follows: 

V£ e R, sign(0 = 



+1 if C > 
— 1 otherwise. 



4.3 Proximity operator of the temporal regularization function 

A simple expression of the proximity operator of function h is not available. We thus propose to 
split this regularization term as a sum of two more tractable functions hi and h^ 

N r /2 

hi(o = kJ2 n r *c 2t - T *c 2t_1 iip (i6) 
t=i 

and 

N r /2-l 

MC) = « E ||T*C 2t+1 -T*C 2 *||^- (17) 

t=l 

Since /ii (respectively /i 2 ) is separable w.r.t the time variable £, its proximity operator can easily 
be calculated based on the proximity operator of each of the involved terms in the sum of Eq. (16) 
(respectively Eq. (17)). 

Indeed, let us consider the following function 

$:C x xC K ^R (18) 

(C*, C'" 1 ) ^ «||t*c' - t^C'X = V o H((\ C'" 1 ), 

where ip(-) = ft||T* • \\p and if is the linear operator defined as 

H : C K x C x — > C K (19) 
(a, b) a — b. 

Its associated adjoint operator if* is therefore given by 

H * : C K — > C K x C x (20) 
a i-)- (a, —a). 

Since we have ifif* = 2Id, the proximity operator of \I/ can easily be calculated using [50, 
Prop. 11]: 

prox^ = prox^ oif = Id + -if* o (prox 2V , - Id) o if. (21) 
The calculation of prox 2 ^ is discussed in [47]. 
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4.4 Parallel Proximal Algorithm (PPXA) 

The function to be minimized has been reexpressed as 

N r 

^st(C) =E E HdV) - S(r)(T*C t )(r)\\l. 1 +g(C t ) + h 1 (C t ) + h 2 (C t ). 

t=l r6{l,...,X}x{l,...,y/fl}x{l,...,Z} 

(22) 

Since Jst is made up of more than two non-necessarily differentiable terms, an appropriate 
solution for minimizing such an optimality criterion is PPXA [51]. In particular, it is important 
to note that this algorithm does not require subiterations as was the case for the constrained 
optimization algorithm proposed in [29]. In addition, the computations in this algorithm can be 
performed in a parallel manner and the convergence of the algorithm to an optimal solution to the 
minimization problem is guaranteed. 

The resulting algorithm for the minimization of the optimality criterion in Eq. (22) is given in 
Algorithm 1. In this algorithm, the weights uji have been fixed to 1/4 for every i E {1, . . . , 4}. The 
parameter 7 has been set to 200 since this value was observed to lead to the fastest convergence in 
practice. The stopping parameter e has been set to 10 -4 . Using these parameters, the algorithm 
typically converges in less than 50 iterations. 

5 Experimental validation in neuroimaging 

This section is dedicated to the experimental validation of the reconstruction algorithm we proposed 
in Section 3.2. Experiments have been conducted on both anatomical and functional data in 
Sections 5.1 and 5.2, respectively. For anatomical data, the proposed 3D-UWR-SENSE algorithm 
(4D-UWR-SENSE without temporal regularization) is compared to the Siemens reconstruction 
pipeline. As regards fMRI validation, results of subject and group-level fMRI statistical analyses 
are compared for two reconstruction pipelines: the one available on the Siemens workstation and 
our own pipeline involving for the sake of completeness either the early UWR-SENSE [29] or the 
4D-UWR-SENSE version of the proposed pMRI reconstruction algorithm. 

5.1 Anatomical images 

For validation purpose, we acquired anatomical and fMRI data on a 3T Siemens Trio magnet. 
Anatomical data has been acquired using a 3D Tl-weighted MP-RAGE pulse sequence at a 1 x 1 x 
1.1 mm 3 spatial resolution (TE — 2.98 ms, TR — 2300 ms, slice thickness = 1.1 mm, transversal 
orientation, FOV = 256 x 240 x 176 mm 3 ). 

To compare the proposed approach to the mSENSE 2 one, Fig. 1 illustrates coronal and axial 
anatomical slices reconstructed with both algorithms while turning off the temporal regularization 
in 4D-UWR-SENSE, so resulting in the so-called 3D-UWR-SENSE approach. Red ellipsoids clearly 
show reconstruction artifacts in the mSENSE reconstruction, which have been removed using our 3D- 
UWR-SENSE approach. It is also worth noticing that similar results have been obtained of different 
subjects. 

In this context, comparisons with other methods (e.g. H-FOCUSS-like methods) are not possible 
2 SENSE reconstruction implemented by the Siemens scanner 
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Algorithm 1 4D-UWR-SENSE: spatio-temporal regularized reconstruction. 

Set ( 7) e) e]0,+oo[ 2 , (wi)i<i<4 €]0,1[ 4 such that Et=i"i = 1, (C| n) )i<i<4 e (C^^f where c\ n) = 
(C^ n) ,C* An \...,C? rAn) ), n = 0, and = ((<;.<»>), ((C^oeo.i^J for every < G {!,... ,4} and 

t £ {1, . . . , AT r }. Set also C (n) = Ei=i f>iCj n) and j( n ) = 0. 
1: repeat 

2: S etpi' (n) =C 4 1,(n) - 
3: for t = 1 to 7V r do 

4: Compute p*' (n) = prox 7 ^ wLg/u;i (Cj' (n) ). 

5: Compute p 2 ' (n) = (prox 7$a/u;2 (C^), (prox^^ . /w2 (C^oj))oeO,i<j<j max ) • 

6: if t is even then 

7: Compute (Ps (n \pl~ 1,(n) ) = Wox^^^/^, Cl~ hM ) 

8: else if t is odd and t > 1 then 

9: Compute (n) ,^ _1 ' (n) ) = prox^/^ (C4 (n) > d" 1 '^)- 

10: end if 
11: if t > 1 then 

12: Set^- 1 'W=Eti^" 1,(n) - 

13: end if 

14: end for 

15: Setpf' (n) =cf' (n) . 

16: Compute P^-W = £? =1 Wi pf- W . 

17: Set = (P 1 '^), W n ),...,P^'( n )). 

18: Set A n e [0,2]. 
19: for z = 1 to 4 do 

20: Set^ = (^^\^^\...,pf-^). 

21: Compute = + A n (2P(") - C (n) -p, (n) ). 

22: end for 

23: Compute C (n+1) = C (n) + A n (P (n) - C (n) ). 

24: ra<-rc+l. 

25: until |Jsx(C (n) ) - JsT(C (n_1) )l < eJsT(C (n_1) ). 

26: Set C = C (n) - 

27: return p £ = T*C £ for every £ G {1, . . . , iV r }. 



since those are only able to deal with 2D+t data. This shows the versatility of our approach which 
allows one, in addition to 2D+t data, to process data collected by 3D acquisition sequences. 

5.2 Functional datasets 

For fMRI data, a Gradient-Echo EPI (GE-EPI) sequence has been used (TE = 30 ms, TR = 
2400 ms, slice thickness = 3 mm, transversal orientation, FOV = 192 mm 2 ) during a cognitive 
localizer [52] protocol. 

This experiment has been designed to map auditory, visual and motor brain functions as well 
as higher cognitive tasks such as number processing and language comprehension (listening and 
reading). It consisted of a single session of N r = 128 scans. The paradigm was a fast event-related 
design comprising sixty auditory, visual and motor stimuli, defined in ten experimental conditions 
(auditory and visual sentences, auditory and visual calculations, left/right auditory and visual 
clicks, horizontal and vertical checkerboards). An L = 32 channel coil was used to enable parallel 
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mSENSE 3D-UWR-SENSE 



Coronal 





Axial 





Figure 1: Coronal and Axial reconstructed slices using mSENSE and 3D-UWR-SENSE (4D-UWR- 
SENSE without temporal regularization) for R — 4 with 1 x 1 x 1.1 mm 3 spatial resolution. Red 
ellipsoids indicate the position of reconstruction artifacts using mSENSE. 

imaging. Ethics approval was given by the local research ethics committee, and fifteen subjects 
gave written informed consent for participation. 



5.2.1 FMRI reconstruction pipeline 



For each subject, fMRI data were collected at the 2x2 mm 2 spatial in-plane resolution using 
different reduction factors (R = 2 or R = 4). Based on the raw data files delivered by the scanner, 
reduced FOV EPI images were reconstructed as detailed in Fig. 2. This reconstruction is performed 
in two stages: 

i) the k-space regridding to account for the non-uniform fc-space sampling during readout 
gradient ramp, which occurs in fast MRI sequences like GE-EPI; 

ii) the Nyquist ghosting correction to remove the odd-even echo inconsistencies during /c-space 
acquisition of EPI images. 



Siemens MRI 


Reading raw FID 




Regridding 


scanner 


data 







Nyquist ghosting 




IFFT 


Reduced FOV 




correction 




images 



Figure 2: Reconstruction pipeline of reduced FOV EPI images from the raw FID data. 
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It must be emphasized here that since no interleaved fc-space sampling is performed during 
the acquisition, and since the central lines of the fc-space are not acquired for each TR due to the 
available imaging sequences on the Siemens scanner, H-FOCUSS-like methods are not applicable 
on the available dataset. 

Once the reduced FOV images are available, the proposed pMRI 4D-UWR-SENSE algorithm and its 
early UWR-SENSE version have been utilized in a final step to reconstruct the full FOV EPI images 
and compared to the mSENSE Siemens solution. For the wavelet-based regularization, dyadic (M = 
2) Symmlet orthonormal wavelet bases [53] associated with filters of length 8 have been used over 
jmax = 3 resolution levels. The reconstructed EPI images then enter in our fMRI study in order 
to measure the impact of the reconstructor choice on brain activity detection. Note also that the 
proposed reconstruction algorithm requires the estimation of the coil sensitivity maps (matrix S(-) 
in Eq. (2)). As proposed in [4], the latter were estimated by dividing the coil-specific images by 
the module of the Sum Of Squares (SOS) images, which are computed from the specific acquisition 
of the fc-space centre (24 lines) before the N r scans. The same sensitivity map estimation is 
then used for all compared methods. Fig. 3 compares the two pMRI reconstruction algorithms to 
illustrate on axial, coronal and sagittal EPI slices how the mSENSE reconstruction artifacts have been 
removed using the 4D-UWR-SENSE approach. The mSENSE reconstructed images acutally present 
large artifacts located both at the centre and boundaries of the brain in sensory and cognitive 
regions (temporal lobes, frontal and motor cortices, ...). This results in SNR loss and thus may 
have a dramatic impact for activation detection in these brain regions. Note that these conclusions 
are reproducible across subjects although the artifacts may appear on different slices (see red circles 
in Fig. 3). It is also worth mentioning that in contrast with the Siemens reconstructor our pipeline 
does not involve any spatial filtering step to improve signal homogeneity across the brain: this 
explains why the images shown in Fig. 3 and delivered by our 4D-UWR-SENSE algorithm seem 
less homogeneous. However, bias field correction can be applied if necessary using specific tools 
such as those available in BrainVISA 3 . 

Regarding computational load, the mSENSE algorithm is carried out on-line and remains 
compatible with real time processing. On the other hand, our pipeline is carried out off-line and 
requires more computations. For illustration purpose, on a biprocessor quadcore Intel Xeon 
CPU® 2.67GHz, one EPI slice is reconstructed in 4 s using the UWR-SENSE algorithm. Using par- 
allel computing strategy and multithreading (through the OMP library) , each EPI volume consisting 
of 40 slices is reconstruced in 22 s. This makes the whole series of 128 EPI images available in about 
47 min. By contrast, the proposed 4D-UWR-SENSE achieves the reconstruction of the series in 
about 40 min, but requires larger memory space due to large data volume processed simultaneously. 



5.2.2 fMRI data pre- processings 

Irrespective of the reconstruction pipeline, the full FOV fMRI images were then prepro- 
cessed using the SPM5 software 4 : preprocessing involves realignment, correction for motion 
and differences in slice acquisition time, spatial normalization, and smoothing with an isotropic 

3 http://brainvisa.info. 

4 htt p : / / www . fil . ion . ucl .ac.uk/spm/ software / spm5 
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Figure 3: Axial, Coronal and Sagittal reconstructed slices using mSENSE and 4D-UWR-SENSE 
for R = 2 and R = 4 with 2x2 mm 2 in-plane spatial resolution. Red circles and ellipsoids indicate 
the position of reconstruction artifacts using mSENSE. 



Gaussian kernel of 4mm full- width at half-maximum. Anatomical normalization to MNI space was 
performed by coregistration of the functional images with the anatomical Tl scan acquired with 
the thirty-two channels head coil. Parameters for the normalization to MNI space were estimated 
by normalizing this scan to the Tl MNI template provided by SPM5, and were subsequently 
applied to all functional images. 



5.2.3 Subject-level analysis 



A General Linear Model (GLM) was constructed to capture stimulus-related BOLD response. 
As shown in Fig. 4, the design matrix relies on ten experimental conditions and thus made up of 
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twenty one regressors corresponding to stick functions convolved with the canonical Haemodynamic 
Response Function (HRF) and its first temporal derivative, the last regressor modelling the baseline. 
This GLM was then fitted to the same acquired images but reconstructed using either the Siemens 
reconstructor or our own pipeline, which in the following is derived from the early UWR-SENSE 
method [29] and from its 4D-UWR-SENSE extension we propose here. Here, contrast estimated 




Figure 4: (a): Design matrix and the Lc-Rc contrast involving two conditions (grouping auditory 
and visual modalities); (b): design matrix and the aC-aS contrast involving four conditions (sen- 
tence, computation, left click, right click). 

images for motor responses and higher cognitive functions (computation, language) were subjected 
to further analyses at the subject and group levels. These two analyses are complementary since 
the expected activations lie in different brain regions and thus can be differentially corrupted by 
reconstruction artifacts as outlined in Fig. 3. More precisely, we studied: 

• the Auditory computation vs. Auditory sentence (aC-aS) contrast which is supposed 
to elicit evoked activity in the frontal and parietal lobes, since solving mental arithmetic task 
involves working memory and more specifically the intra-parietal sulcus [54]: see Fig. 4(b); 

• the Left click vs. Right click (Lc-Rc) contrast for which we expect evoked activity in 
the right motor cortex (precentral gyrus, middle frontal gyrus). Indeed, the Lc-Rc contrast 
defines a compound comparison involving two motor stimuli which are presented either in the 
visual or auditory modality. This comparison aims therefore at detecting lateralization effect 
in the motor cortex: see Fig. 4(a). 

Interestingly, these two contrasts were chosen because they summarized well different situa- 
tions (large vs small activation clusters, distributed vs focal activation pattern, bilateral vs uni- 
lateral activity) that occurred for this paradigm when looking at sensory areas (visual, auditory, 
motor) or regions involved in higher cognitive functions (reading, calculation). In the following, our 
results are reported in terms of Student-t maps thresholded at a cluster-level p = 0.05 corrected for 
multiple comparisons according to the Family Wise Error Rate (FWER) [55,56]. Complementary 
statistical tables provide corrected cluster and voxel-level p-values, maximal t-scores and corre- 
sponding peak positions both for R = 2 and R — 4. Note that clusters are listed in a decreasing 
order of significance. 
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Concerning the aC-aS contrast, Fig. 5 [top] shows for the most significant slice and R — 2 that 
all pMRI reconstruction algorithms succeed in finding evoked activity in the left parietal and frontal 
cortices, more precisely in the inferior parietal lobule and middle frontal gyrus according to the AAL 
template 5 . Table 1 also confirms a bilateral activity pattern in parietal regions for R — 2. Moreover, 
for R = 4, Fig. 5 [bottom] illustrates that our pipeline (UWR-SENSE and 4D-UWR-SENSE) and 
preferentially the proposed 4D-UWR-SENSE scheme enables to retrieve reliable frontal activity 
elicited by mental calculation, which is lost by the the mSENSE algorithm. From a quantitative 
viewpoint, the proposed 4D-UWR-SENSE algorithm finds larger clusters whose local maxima are 
more significant than the ones obtained using mSENSE and UWR-SENSE, as reported in Table 1. 
Concerning the most significant cluster for R = 2, the peak positions remain stable whatever the 
reconstruction algorithm. However, examining their significance level, one can first measure the 
benefits of wavelet-based regularization when comparing UWR-SENSE with mSENSE results and 
then additional positive effects of temporal regularization and 3D wavelet decomposition when 
looking at the 4D-UWR-SENSE results. These benefits are also demonstrated for R = 4. 
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Figure 5: Subject-level student-t maps superimposed to anatomical MRI for the aC-aS contrast. 
Data have been reconstructed using the mSENSE, UWR-SENSE and 4D-UWR-SENSE, respectively. 
Neurological convention: left is left. The blue cross shows the maximum activation peak. 

Fig. 6 illustrates another property of the proposed pMRI pipeline, i.e. its robustness to the 
between-subject variability. Indeed, when comparing subject-level student-t maps reconstructed 
using the different pipelines (R = 2), it can be observed that the mSENSE algorithm fails to detect 
any activation cluster in the expected regions for the second subject (see Fig. 6 [bottom]). By 
contrast, our 4D-UWR-SENSE method retrieves more coherent activity while not exactly at the 
same position as for the first subject. 

For the Lc-Rc contrast on the data acquired with i? = 2, Fig. 7 [top] shows that all reconstruction 
methods enable to retrieve the expected activation in the right precentral gyrus. However, when 



3 available in the xjView toolbox of SPM5. 
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Table 1: Significant statistical results at the subject-level for the aC-aS contrast (corrected for 
multiple comparisons at p = 0.05). Images were reconstructed using the mSENSE, UWR-SENSE 
and 4D-UWR-SENSE algorithm for R = 2 and R = 4. 
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voxel-level 
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R = 2 


mSENSE 


< 10~ 3 


320 


< 10~ 3 


6.40 


-32 -76 45 


< 10~ 3 


163 


< 10~ 3 


5.96 


-4 -70 54 


< 10~ 3 


121 


< 10~ 3 


6.34 


34 -74 39 


< 10~ 3 


94 


< 10~ 3 


6.83 


-38 4 24 


UWR-SENSE 


< 10~ 3 


407 


< 10~ 3 


6.59 


-32 -76 45 


< 10~ 3 


164 


< 10~ 3 


5.69 


-6 -70 54 


< 10" 3 


159 


< 10" 3 


5.84 


32 -70 39 


< 10" 3 


155 


< 10" 3 


6.87 


-44 4 24 


4D-UWR-SENSE 


< 10" 3 


454 


< 10" 3 


6.54 


-32 -76 45 


< 10" 3 


199 


< 10" 3 


5.43 


-6 26 21 


< 10" 3 


183 


< 10" 3 


5.89 


32 -70 39 


< 10" 3 


170 


< 10" 3 


6.90 


-44 4 24 


R = 4 


mSENSE 


< 10" 3 


58 


0.028 


5.16 


-30 -72 48 


4D-UWR-SENSE 


< 10" 3 


94 


0.003 


5.91 


-32 -70 48 


< icr 3 


60 


0.044 


4.42 


-6 -72 54 


4D-UWR-SENSE 


< 10" 3 


152 


< 10" 3 


6.36 


-32 -70 48 


< icr 3 


36 


0.009 


5.01 


-4 -78 48 


< 10" 3 


29 


0.004 


5.30 


-34 6 27 



looking more carefully at the statistical results (see Table 2), our pipeline and more preferentially 
the 4D-UWR-SENSE algorithm retrieves an additional cluster in the right middle frontal gyrus. 
On data acquired with R = 4, the same Lc-Rc contrast elicits similar activations, i.e. in the same 
region. As demonstrated in Fig. 7[bottom], this activity is enhanced when pMRI reconstruction 
is performed with our pipeline. Quantitative results in Table 2 confirm numerically what can be 
observed in Fig. 7: larger clusters with higher local t-scores are detected using the 4D-UWR-SENSE 
algorithm, both for R = 2 and R = 4. Also, a larger number of clusters is retrieved for R = 2 using 
wavelet-based regularization. 

Fig. 8 reports on the robustness of the proposed pMRI pipeline to the between-subject vari- 
ability for this motor contrast. Since sensory functions are expected to generate larger BOLD 
effects (higher SNR) and appear more stable, our comparison takes place at R = 4. Two subject- 
level student-t maps reconstructed using the different pMRI algorithms are compared in Fig. 8. For 
the second subject, one can observe that the mSENSE algorithm fails to detect any activation clus- 
ter in the right motor cortex. By contrast, our 4D-UWR-SENSE method retrieves more coherent 
activity for this second subject in the expected region. 
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Figure 6: Between-subject variability of detected activation for the aC-aS contrast at R 
Neurological convention. The blue cross shows the activation peak. 
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Figure 7: Subject-level student-t maps superimposed to anatomical MRI for the Lc-Rc contrast. 
Data have been reconstructed using the mSENSE, UWR-SENSE and 4D-UWR-SENSE, respectively. 
Neurological convention. The blue cross shows the activation peak. 



To summarize, for these two contrasts our 4D-UWR-SENSE algorithm always outperforms the 
alternative reconstruction methods in terms of statistical significance (number of clusters, cluster 
extent, peak values,...) but also in terms of robustness. 
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Table 2: Significant statistical results at the subject-level for the Lc-Rc contrast (corrected for 
multiple comparisons at p = 0.05). Images were reconstructed using the mSENSE, UWR-SENSE 
and 4D-UWR-SENSE algorithms for R = 2 and R = 4. 





cluster-level 


voxel-level 


p- value 


Size 


p- value 


T-score 
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R = 2 


mSENSE 


< 10" 3 


79 


< 10" 3 


6.49 


38 -26 66 


UWR-SENSE 


< 10" 3 


144 


0.004 


5.82 


40 -22 63 


0.03 


21 


0.064 


4.19 


24 -8 63 


4D-UWR-SENSE 


< 10" 3 


172 


0.001 


6.78 


34 -24 69 


< 10" 3 


79 


0.001 


6.49 


38 -26 66 


R = 4 


mSENSE 


0.006 
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4.82 
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< 10" 3 
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5.57 


40 -24 66 
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Figure 8: Between-subject variability of detected activation for the Lc-Rc contrast at R = 4. 
Neurological convention. The blue cross shows the activation peak. 



5.2.4 Group-level analysis 



Due to between-subject anatomical and functional variability, group-level analysis is necessary 
in order to derive robust and reproducible conclusions at the population level. For this validation, 
random effect analyses (RFX) involving fifteen healthy subjects have been conducted on the contrast 
maps we previously investigated at the subject level. More precisely, one-sample Student-t test was 
performed on the subject-level contrast images (eg, Lc-Rc, aC-aS,... images) using SPM5. 

For the aC-aS contrast, Maximum Intensity Projection (MIP) student-t maps are shown in 
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Figure 9: Group-level student-t maps for the aC-aS contrast where data have been reconstructed 
using the mSENSE, UWR-SENSE and 4D-UWR-SENSE for R = 2 and R = 4. Neurological conven- 
tion. Red arrows indicate the global maximum activation peak. 



Fig. 9. First, they illustrate that irrespective of the reconstruction method larger and more signifi- 
cant activations are found on datasets acquired with R — 2 providing the better SNR. Second, for 
R — 2, visual inspection of Fig. 9[top] confirms that only the 4D-UWR-SENSE algorithm allows us 
to retrieve significant bilateral activations in the parietal cortices (see axial MIP slices) in addition 
to larger cluster extent and a gain in significance level for the stable clusters across the differ- 
ent reconstruct ors. Similar conclusions can be drawn when looking at Fig. 9 [bottom] for R — 4. 
Complementary results are available in Table 3 for R — 2 and R — 4 and allow us to numerically 
validate this visual comparison: 

• Whatever the reconstruction method in use, the statistical performance is much more signif- 
icant using R — 2, especially at the cluster level since the cluster extent decreases by one 
order of magnitude. 

• Voxel and cluster-level results are enhanced using the 4D-UWR-SENSE approach instead of 
the mSENSE reconstruction or its early UWR-SENSE version. 

Fig. 10 reports similar group-level MIP results for R — 2 and R — 4 concerning the Lc-Rc 
contrast. It is shown that whatever the acceleration factor R in use, our pipeline enables to detect 
a much more spatially extended activation area in the motor cortex. This visual inspection is quan- 
titatively confirmed in Table 4 when comparing the detected clusters using our 4D-UWR-SENSE 
approach with those found by mSENSE, again irrespective of R. Finally, the 4D-UWR-SENSE 
algorithm outperforms the UWR-SENSE one, which corroborates the benefits of the proposed 
spatio-temporal regularization scheme. 
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Table 3: Significant statistical results at the group-level for the aC-aS contrast (corrected for 
multiple comparisons at p = 0.05). Images were reconstructed using the mSENSE, UWR-SENSE 
and 4D-UWR-SENSE algorithms for R = 2 and R = 4. 





cluster-level 


voxel-level 


p- value 


Size 


p- value 


T-score 


Position 


R = 2 


mSENSE 


< 10~ 3 


361 


0.014 


7.68 


-6 -22 45 


< 10~ 3 


331 


0.014 


8.23 


-40 -38 42 


< 10~ 3 


70 


0.014 


7.84 


-44 6 27 


UWR-SENSE 


< icr 3 


361 


0.014 


7.68 


-6 22 45 


< icr 3 


331 


0.014 


7.68 


-44 -38 42 


< icr 3 


70 


0.014 


7.84 


-44 6 27 


4D-UWR-SENSE 


< icr 3 


441 


< 10" 3 


9.45 


-32 -50 45 


< icr 3 


338 


< 10" 3 


9.37 


-6 12 45 


< 10" 3 


152 


0.010 


7.19 


30 -64 48 


R = 4 


mSENSE 


0.003 


14 


0.737 


5.13 


-38 -42 51 


UWR-SENSE 


< 10" 3 


41 


0.274 


5.78 


-50 -38 -48 


< 10" 3 


32 


0.274 


5.91 


2 12 54 


4D-UWR-SENSE 


< 10" 3 


37 


0.268 


6.46 


-40 -40 54 


< 10" 3 


25 


0.268 


6.37 


-38 -42 36 


< 10" 3 
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Figure 10: Group-level student-t maps for the Lc-Rc contrast where data have been reconstructed 
using the mSENSE, UWR-SENSE and 4D-UWR-SENSE for R = 2 and R = 4. Neurological conven- 
tion. Red arrows indicate the global maximum activation peak. 



6 Discussion and conclusion 

The contribution of the present paper is twofold. First, we proposed a novel reconstruction method 
that relies on a 3D wavelet transform and accounts for temporal dependencies in successive fMRI 
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Table 4: Significant statistical results at the group-level for the Lc-Rc contrast (corrected for 
multiple comparisons at p = 0.05). Images were reconstructed using the mSENSE, UWR-SENSE 
and 4D-UWR-SENSE algorithms for R = 2 and R = 4. 





cluster-level 


voxel-level 


p- value 


Size 


p-value 


T-score 


Position 


R = 2 


mSENSE 


< 10" 3 


354 


< 10" 3 


9.48 


38 -22 54 


0.001 


44 


0.665 


6.09 


-4 -68 -24 


UWR-SENSE 


< 10" 3 


350 


0.005 


9.83 


36 -22 57 


< 10" 3 


35 


0.286 


7.02 


4 -12 51 


4D-UWR-SENSE 


< 10" 3 


377 


0.001 


11.34 


36 -22 57 


< 10" 3 


53 


< 10" 3 


7.50 


8 -14 51 


< 10" 3 


47 


< 10" 3 


7.24 


-18 -54 -18 


R = A 


mSENSE 


< 10" 3 


38 


0.990 


5.97 


32 -20 45 


UWR-SENSE 


< 10" 3 


163 


0.128 


7.51 


46 -18 60 


4D-UWR-SENSE 


< 10" 3 


180 


0.111 


7.61 


46 -18 60 



volumes. As a particular case, the proposed method allows us to deal with 3D acquired anatom- 
ical data when no temporal repetition is considered. Second, when artifacts are superimposed to 
brain activation, we showed that the choice of the parallel imaging reconstruction algorithm has a 
significant influence on the statistical sensitivity in fMRI at the subject and group-levels and may 
enable whole brain neuroscience studies at high spatial resolution. 

Practically speaking, we showed that whole brain acquisition can be routinely used at a spatial 
in-plane resolution of 2 x 2mm 2 in a short and constant repetition time (TR = 2.4s) provided that 
a reliable pMRI reconstruction pipeline is chosen. In this paper, we demonstrated that our 4D- 
UWR-SENSE reconstruction algorithm meets this goal. To draw this conclusion, our comparison 
took place at the statistical analysis level and relied on quantitative criteria (voxel- and cluster-level 
corrected p-values, t-scores, peak positions) both at the subject and group levels. In particular, 
we showed that our 4D-UWR-SENSE approach outperforms both its UWR-SENSE ancestor [29] 
and the Siemens mSENSE reconstruction in terms of statistical significance and robustness. This 
emphasized the benefits of combining temporal and 3D regularizations in the wavelet domain. The 
usefulness of 3D regularization in reconstructing 3D anatomical images was also shown, especially 
in more degraded situations (R = 4) where regularization plays a prominent role. The validity of 
our conclusions lies in the reasonable size of our datasets since the same cohort participated to 
parallel imaging acquisitions using two different acceleration factors (R = 2 and R = 4). 

At the considered spatio-temporal compromise (2x2x3mm 3 and TR = 2.4 s), we also illustrated 
the impact of increasing the acceleration factor (passing from R = 2 to R = 4) on the statistical 
sensitivity at the subject and group levels for a given reconstruction algorithm. We performed 
this comparison to anticipate what could be the statistical performance for detecting evoked brain 
activity on data requiring this acceleration factor, such as high spatial resolution EPI images (eg 
1.5 x 1.5mm 2 in-plane resolution ones) acquired at the same short TR. Our conclusions were 
balanced depending on the contrast of interest: when looking at the aC-aS contrast, involved in 
the fronto-parietal circuit, it turned out that R = 4 was not reliable enough to recover significant 
group-level activity at 3 Tesla: the SNR loss was too important and should be compensated by 
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an increase of the static magnetic field (e.g. passing from 3 to 7 Tesla). However, the situation 
becomes acceptable for the Lc-Rc motor contrast, which elicits activation in motor regions: our 
results brought evidence that the 4D-UWR-SENSE approach enables the use of R — 4 for this 
contrast. 

To summarize, the compromise between the acceleration factor and spatial in-plane resolution 
should be selected with care depending on the regions involved in the fMRI paradigm. As a 
consequence, high resolution fMRI studies can be conducted using high speed acquisition (short 
TR and large R value) provided that the expected BOLD effect is strong, as experienced in primary 
motor, visual and auditory cortices. Of course, the use of an efficient reconstruction method such 
as the one proposed is a pre-requisite to shift this compromise towards larger R values and higher 
spatial resolution and it could be optimally combined with ultra high magnetic fields. 

A direct extension of the present work, which is actually in progress, consists of studying the 
impact of tight frames instead of wavelet bases to define more suitable 3D transforms. However, 
unsupervised reconstruction becomes more challenging in this framework since the estimation 
of hyper-parameters becomes cumbersome (see [57] for details). Ongoing work will concern the 
combination of the present contribution with the joint detection estimation approach of evoked 
activity [58, 59] to go beyond the GLM framework and to evaluate how the pMRI reconstruction 
algorithm also impacts HRF estimation. Another extension of our work would concern the 
combination of our wavelet-regularized reconstruction with the WSPM approach [60] in which 
statistical analysis is directly performed in the wavelet transform domain. 

Appendix 

A Maximum likelihood estimation of regularization parameters 

A rigorous way of addressing the regularization parameter choice would be to consider that the sum 
of the regularization functions g and h corresponds to the minus-log-likelihood of a prior distribution 

/(.;©) where = (M ajmax , a ajmax , /3 aJmax , (t*oj,aoj,P j) oe oi<j<j max > K >P)> and to maximize 
the integrated likelihood of the data. This would however entail two main difficulties. On the one 
hand, this would require to integrate out the sought image decomposition ( and to iterate between 
image reconstruction and hyper-parameter estimation. Methods allowing us to perform this task 
are computationally intensive [61]. On the second hand, the partition function of the distribution 
/(•; 0) does not take a closed form and we would thus need to resort to numerical methods [62-64] 
to compute it. To alleviate the computational burden, akin to [65] we shall proceed differently by 
assuming that a reference full FOV image p is available, and so is its wavelet decomposition £ = Tp. 
In practice, our reference image p is obtained using ID-SENSE reconstruction at the same R value. 
We then apply an approximate ML procedure which consists of estimating separately the spatial 
and temporal parameters. Although this approach is not optimal from a theoretical standpoint, it 
is quite simple and it was observed to provide satisfactory results in practice. Alternative solutions 
based on Monte Carlo methods [57] or Stein's principle [66] can also be thought of, at the expense 
of an additional computational complexity. 
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A.l Spatial regularization parameters 

For the spatial hyper-parameter estimation task, we will assume that the real and imaginary parts 
of the wavelet coefficients 6 are modelled by the following Generalized Gauss-Laplace (GGL) distri- 
bution: 

V^R, mwft-y/JL __ . (23) 

For each resolution level j and orientation o, a^j and are estimated from C i0 - as follows 
(we proceed similarly to estimate /I*™, S*™ and /3^- by replacing Re(-) by Im(-)): 

(jug, ^fjjfj) = argmax /(Re(C oJ ); a, /3) 

0,a,/3)£]RxM + xR^ 

= argmax V log /(Re(Coj,ife); a, P) 

0,a,/3)eRxR+xR^ k=1 

« |Re(C ,i,fc - n)\+- V |Re(Co,i,fc - A*)| 2 

^ )U ,/ J ;CiB.AJ B . + AiB. + fe=1 Z fe=1 

+ - if log/3 + Kj log ( erfc( ^ } ) }• (24) 

This three-dimensional minimization problem does not admit a closed form solution. Hence, we 
can compute the ML estimated parameters using the zero-order Powell optimization method [67] . 



A. 2 Temporal regularization parameter 

For the temporal hyper-parameter estimation task, we will assume that, at a given voxel, the 
temporal noise is distributed according to the following generalized Gaussian (GG) distribution: 

VeER, /(q^P) = P 2r(1/p) • (25) 

Akin to the spatial hyper-parameter estimation, reference images {^)i<t<N r are made available 
based on a ID-SENSE reconstruction, where Vt E {l,...,7V r }, = T*£*. We consider that 
at spatial position r, the temporal noise vector e r = [/^(r) — p 1 (r),p 3 (r) — /^(r), . . . ,p Nr (r) — 
p- /Vr_1 (r)] T is a realization of a full independent GG prior distribution and we adjust the temporal 
hyper-parameter vector directly from it. It should be noted here that the considered model 

for the temporal noise accounts for correlations between successive observations usually considered 
in the fMRI literature. It also presents more flexibility than the Gaussian model, which corresponds 



6 A similar approach is adopted for the approximation coefficients. 



25 



to the particular case when p = 2. Estimates k and p of the parameters are then obtained as follows: 
(n,p) = argmax f{e r \K,p) 

(k,p)GM+x[1,+cx)[ 

= argmax log / (e r ; K,p) 

(«,p)GM+x[l,+cx)[ 



= argmin « |^(r) - ?(r)|* - (tf r - 1) log (-^— ) . (26) 



Note that in the above minimization, for a given value of p, the optimal value of k admits the 
following closed form: 

8= „ 1 JV --' . (27) 

A zero-order Powell optimization method can then be used to solve the resulting one-variable 
minimization problem. To reduce the computational complexity of this estimation, it is only 
performed on the brain mask, and the temporal regularization parameter n is set to zero for voxels 
belonging to the image background. 
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